
# Standaard libaries voor visualiseren en dataverwerking
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns

# Libaries gebruikt bij de seizoens analyse
from statsmodels.tsa.seasonal import seasonal_decompose
from scipy.signal import find_peaks
from statsmodels.tsa.deterministic import CalendarFourier
from statsmodels.tsa.deterministic import DeterministicProcess

# Libaries gebruikt bij trend analyse

# Libaries gebruikt bij autocorrelatie
import statsmodels.api as sm
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf

# Libaries gebruikt bij functies rondom modelleren
import datetime
from sklearn.metrics import mean_squared_error
from sklearn.model_selection import GridSearchCV
from sklearn.model_selection import train_test_split
from sklearn.model_selection import TimeSeriesSplit
from statsmodels.tsa.stattools import adfuller

# Libaries van de modellen
from sklearn.linear_model import LinearRegression
from sklearn.neighbors import KNeighborsRegressor
from sklearn.tree import DecisionTreeRegressor
from sklearn.ensemble import RandomForestRegressor
from sklearn.ensemble import StackingRegressor
from xgboost import XGBRegressor
from statsmodels.tsa.stattools import adfuller
from statsmodels.tsa.statespace.sarimax import SARIMAX
from prophet import Prophet


# Inladen van de test dataset
data_test = pd.read_csv('test.csv', sep = ',')

# Inladen van de train dataset
data_train = pd.read_csv('train.csv', sep = ',')


# Tonen van de eerste 5 regels
data_train.head()


# Tonen van de kern informatie over kolommen
data_train.info()

<class 'pandas.core.frame.DataFrame'>
RangeIndex: 16637 entries, 0 to 16636
Data columns (total 8 columns):
 #   Column      Non-Null Count  Dtype  
---  ------      --------------  -----  
 0   date_hour   16637 non-null  object 
 1   holiday     16637 non-null  int64  
 2   weathersit  16637 non-null  int64  
 3   temp        16637 non-null  float64
 4   atemp       16637 non-null  float64
 5   hum         16637 non-null  float64
 6   windspeed   16637 non-null  float64
 7   cnt         16637 non-null  int64  
dtypes: float64(4), int64(3), object(1)
memory usage: 1.0+ MB


data_train['date_hour'] = pd.to_datetime(data_train['date_hour'])
data_train.info()

<class 'pandas.core.frame.DataFrame'>
RangeIndex: 16637 entries, 0 to 16636
Data columns (total 8 columns):
 #   Column      Non-Null Count  Dtype         
---  ------      --------------  -----         
 0   date_hour   16637 non-null  datetime64[ns]
 1   holiday     16637 non-null  int64         
 2   weathersit  16637 non-null  int64         
 3   temp        16637 non-null  float64       
 4   atemp       16637 non-null  float64       
 5   hum         16637 non-null  float64       
 6   windspeed   16637 non-null  float64       
 7   cnt         16637 non-null  int64         
dtypes: datetime64[ns](1), float64(4), int64(3)
memory usage: 1.0 MB


# Tonen van informatie over numerieke kolommen
data_train.describe()


min_value = data_train['windspeed'].min()
max_value = data_train['windspeed'].max()

data_train['windspeed'] = (data_train['windspeed'] - min_value) / (max_value - min_value)

data_train.describe()


def boxplot_outliers(dataframe, column_name):
    plt.figure(figsize=(8, 6))
    plt.boxplot(dataframe[column_name], vert=False)
    plt.title(column_name)
    plt.show()


for col in data_train.columns:
    if col != 'date_hour':
        boxplot_outliers(data_train, col)


outliers_windspeed = data_train[data_train['windspeed'] > 0.58]

outliers_windspeed


class Grafieken:
    """
    Een class voor eenvoudig visualiseren met Matplotlib
    en seaborn.

    Parameters:
    ----------
    df : pandas.DataFrame
        Een pandas.DataFrame die data bevat om mee te
        visualiseren.
    """

    def __init__(self, df, target=False):
        """
        De constructor voor de Grafieken class

        Parameters:
        ----------
        df : pandas.DataFrame
            Een pandas.DataFrame die data bevat om mee te
            visualiseren.

        target : str (default = False)
            De kolomnaam van de target kolom, indien toepasbaar
        """
        # Aanmaken self.df
        self.df = df.drop('date_hour', axis=1)

        # Zetten van de target
        if target:
            self.target = target

        # Aanmaken grafiek standaarden
        self.figsize = (15,5)
        self.pallete = 'deep'


    def boxplots(self):
        """
        Maakt een aantal box plots gelijk aan het aantal kolommen
        in de Dataframe. De box zijn tegen één specifieke kolom
        opgezet, om de verdeling t.o.v. die kolom te tonen.

        Parameters:
        ----------------
        df : Pandas DataFrame
            Een pandas DataFrame met kolommen waarvan je
            de verdeling wilt onderzoeken.

        kolom : str
            De naam van de kolom waarbij de verdeling bekeken
            word, ten opzichte van de andere kolommen.

        Returns:
        ---------------
        None :
            In plaats van iets te returnen laat de
            functie meerdere boxplots zien.
        """
        # For-loop om elke kolom een eigen boxplot te geven
        for col in self.df.columns:
            # if-statement om te voorkomen dat de
            # kolom-parameter ook een boxplot krijgt
            if col != self.target:
                # Zetten van de grootte van de plot
                plt.subplots(figsize=self.figsize)

                # Maken en benoemen van de assen van de boxplots
                sns.boxplot(x=col, y=self.target,
                            data=self.df, color='pink')
                plt.xlabel(col)
                plt.ylabel(self.target)
                plt.title('{} vs. {}'.format(col, self.target))
                plt.show()
            
    def barplots(self):
        """
        Maakt een aantal bar plots gelijk aan het aantal kolommen
        in de Dataframe. De bar is tegen één specifieke kolom
        opgezet, om de verdeling t.o.v. die kolom te tonen.

        Parameters:
        ----------------
        df : Pandas DataFrame
            Een pandas DataFrame met kolommen waarvan
            je de verdeling wilt onderzoeken.

        kolom : str
            De naam van de kolom waarbij de verdeling
            bekeken word, ten opzichte van de andere kolommen.

        Returns:
        ---------------
        None :
            In plaats van iets te returnen laat
            de functie meerdere barplots zien.
        """
        # For-loop om elke kolom een eigen barplot te geven
        for col in self.df.columns:
            # if-statement om te voorkomen dat de
            # kolom-parameter ook een barplot krijgt
            if col != self.target:
                # Zetten van de grootte van de plot
                plt.subplots(figsize=self.figsize)

                # Maken en benoemen van de assen van de barplot
                sns.barplot(x=col, y=self.target,
                            data=self.df, color= 'skyblue')
                plt.xlabel(col)
                plt.ylabel(self.target)
                plt.title('{} vs. {}'.format(col, self.target))
                plt.show()
                
    def lineplots(self):
        """
        Maakt een aantal line plots gelijk aan het aantal kolommen
        in de Dataframe. De line is tegen één specifieke kolom
        opgezet, om de verdeling t.o.v. die kolom te tonen.

        Parameters:
        ----------------
        df : Pandas DataFrame
            Een pandas DataFrame met kolommen waarvan
            je de verdeling wilt onderzoeken.

        kolom : str
            De naam van de kolom waarbij de verdeling
            bekeken word, ten opzichte van de andere kolommen.

        Returns:
        ---------------
        None :
            In plaats van iets te returnen laat
            de functie meerdere lineplots zien.
        """
        # For-loop om elke kolom een eigen barplot te geven
        for col in self.df.columns:
            # if-statement om te voorkomen dat de
            # kolom-parameter ook een barplot krijgt
            if col != self.target:
                # Zetten van de grootte van de plot
                plt.subplots(figsize=self.figsize)

                # Maken en benoemen van de assen van de barplot
                sns.lineplot(x=col, y=self.target,
                            data=self.df, color= 'green')
                plt.xlabel(col)
                plt.ylabel(self.target)
                plt.title('{} vs. {}'.format(col, self.target))
                plt.show()

    def heatmap(self):
        """
        Deze functie maakt een heatmap op basis
        van de DataFrame waarmee de class is gemaakt
        """
        # Aanmaken Heatmap
        fig, ax = plt.subplots(figsize=(5,15))
        ax = sns.heatmap(self.df.corr()[[self.target]].sort_values(
            by=[self.target], ascending=False
            ),square=False, annot=True, cmap='RdYlGn')
        plt.show()

    def timeplot(self, x_kolom, hue=None):
        """
        Maakt een bar en een boxplot op basis van een gegeven tijd.

        Parameters:
        ----------
        x_kolom : str
            De naam van de tijd gebaseerde kolom

        hue : str
            De naam van de kolom waarop de bars en boxes
            worden ingedeeld

        Returns:
        ----------
        None :
            In plaats van iets te returnen laat
            de functie een bar- en een boxplot zien.
        """
        # Aanmaken van plotgrootte en barplot
        plt.subplots(figsize=self.figsize)
        sns.barplot(data=self.df, x=x_kolom, y=self.target,
                    hue=hue, errorbar=None,
                    palette=self.pallete)

        # Aanmaken van de assen
        plt.title(f'Verhuuraantallen per {x_kolom}')
        plt.ylabel('Verhuuraantallen')
        plt.show()

        # Aanmaken van plotgrootte en boxplot
        plt.subplots(figsize=self.figsize)
        sns.boxplot(data=self.df, x=x_kolom, y=self.target,
                    hue=hue, palette=self.pallete)

        # Aanmaken van de assen
        plt.title(f'Verdeling van Verhuuraantallen per {x_kolom}')
        plt.ylabel('Verhuuraantallen')
        plt.show()


graf = Grafieken(data_train, 'cnt')


graf.heatmap()


graf.boxplots()


graf.barplots()


graf.lineplots()


# Aanmaken van nieuw DataFrame
df_time = data_train[['date_hour', 'cnt']].copy()
df_time['date_hour'] = pd.to_datetime(df_time['date_hour'])
df_time = df_time.set_index('date_hour')
df_time.head()


# Aanmaken van standaard plot
ax = df_time.plot(figsize=(15, 5))

# Toevoegen van shading om beter te zien waar 2012 begint
ax.axvspan('2012-01-01', '2013-01-01', color='green', alpha=0.15)

# Aanpassen van de plot labels en titel
ax.set_xlabel("Datum")
ax.set_ylabel("Verhuuraantallen")
ax.set_title("Verhuuraantallen per uur")

# Tonen van de plot
plt.show()


# Aanmaken van standaard plot
ax = df_time.plot(figsize=(15, 5))

# Toevoegen van shading om beter te zien welke data April is
ax.axvspan('2012-04-01', '2012-05-01', color='green', alpha=0.15)

# Aanpassen van de plot labels en titel
ax.set_xlim('2012-03-22', '2012-05-08')
ax.set_xlabel("Datum")
ax.set_ylabel("Verhuuraantallen")
ax.set_title("Verhuuraantallen per uur van 22 maart tot 8 mei")

# Tonen van de plot
plt.show()


# Gebruiken van seasonal_decompose om seizoens patronen te visualiseren
result = seasonal_decompose(df_time['cnt'], model='additive', period=24*7)

# Maken van de plot
plt.figure(figsize=(15, 5))
result.seasonal.plot()
plt.title('Seizoens Component')
plt.xlim('2012-04-01', '2012-05-01')
plt.xlabel('Datum')
plt.ylabel('Seizoens-waarde')
plt.show()


# Gebruiken van seasonal_decompose om seizoens patronen te visualiseren
result = seasonal_decompose(df_time['cnt'], model='additive', period=24*7)

# Maken van de plot
plt.figure(figsize=(15, 5))
result.seasonal.plot()
plt.title('Seizoens Component')
plt.xlim('2012-04-01', '2012-04-10')
plt.xlabel('Datum')
plt.ylabel('Seizoens-waarde')
plt.show()


# Berekenen aantal datapunten
N = len(df_time)

# Uitvoeren Fourier Transform
ft = np.fft.fft(df_time['cnt'])

# Bereken de kracht van de frequenties
magnitude = 2.0/N * np.abs(ft[:N//2])

# Zoeken en tonen van de pieken
peaks, _ = find_peaks(magnitude, height=75)
print(f"Peaks found at {peaks}")

# Plotten
plt.figure(figsize=(15,5))
plt.plot(magnitude)
plt.xlabel('Frequentie (1/N)')
plt.ylabel('Amplitude')
plt.title('Frequentie domein van Verhuursaantallen')
plt.plot(peaks, magnitude[peaks], ".", markersize=10)

plt.show()

Peaks found at [692 694]


freq = peaks * (1 / N)
periode = 1 / freq
print(f'De periode van FFT is: {np.round(periode, 0)}')

De periode van FFT is: [24. 24.]


# Aanmaken fourier voor W en D
fourier_w = CalendarFourier(freq="W", order=1)
fourier_d = CalendarFourier(freq="D", order=1)

dp = DeterministicProcess(
    index=df_time.index,
    constant=False,
    order=1,
    seasonal=False,
    additional_terms=[fourier_w, fourier_d],
    drop=True,
)

fourier_feat = dp.in_sample()

ax = fourier_feat.plot(y='sin(1,freq=D)',
                       xlim=['2012-01-01', '2012-01-03'])
ax = fourier_feat.plot(ax=ax, y='cos(1,freq=D)',
                       xlim=['2012-01-01', '2012-01-03'])
plt.show()


ax = fourier_feat.plot(y='sin(1,freq=W-SUN)',
                       xlim=('2012-01-01', '2012-01-15'))
ax = fourier_feat.plot(ax=ax, y='cos(1,freq=W-SUN)',
                       xlim=('2012-01-01', '2012-01-15'))
plt.show()


df_trend = data_train.set_index('date_hour')
df_trend


# Index naar datetime
df_trend.index = pd.to_datetime(df_trend.index, format='%Y-%m-%d %H:%M:%S')  
df_trend = df_trend.asfreq('H')  

# Linear regression om de trendlijn te maken
X = pd.to_numeric(df_trend.index).values.reshape(-1, 1)
y = df_trend['cnt'].values
trendmodel = LinearRegression().fit(X, y)

# trendlijn voorspellen
trendline = trendmodel.predict(X)

# Originele timeseries plotten met trendlijn
plt.figure(figsize=(12, 6))
plt.plot(df_trend['cnt'], label='Original timeseries')
plt.plot(df_trend.index, trendline, label='Trendline', color='red')
plt.legend()
plt.title('Timeseries and Trendline')
plt.show()


# These are the "time lags"
shifts = np.arange(1, 25).astype(int)

# Use a dictionary comprehension to create name: value pairs, one pair per shift
shifted_data = {"lag_{}_hour".format(hour_shift): data_train['cnt'].shift(hour_shift) for hour_shift in shifts}

# Convert into a DataFrame for subsequent use
data_train_shifted_data = pd.DataFrame(shifted_data)

# Drop rows with NaN values introduced by the shift operation
data_train_shifted_data.dropna(inplace=True)

# Assuming 'target_variable' is the column you want to predict
target_variable = data_train['cnt'].iloc[data_train_shifted_data.index]

model = LinearRegression() 
model.fit(data_train_shifted_data, target_variable)

fig, ax = plt.subplots()
ax.bar(data_train_shifted_data.columns, model.coef_)
ax.set(xlabel='Coefficient name', ylabel='Coefficient value')
# Set formatting so it looks nice
plt.setp(ax.get_xticklabels(), rotation=45, horizontalalignment='right')

# Plot the first 100 samples of each
ax = data_train_shifted_data.iloc[:100].plot(cmap=plt.cm.viridis)
target_variable.iloc[:100].plot(color='r', lw=2)
ax.legend(loc='best')
plt.show()


fig, (ax1, ax2) = plt.subplots(2,1, figsize=(12,8))
 
# Plot ACF 
plot_acf(data_train['cnt'],lags=24, zero=False, ax=ax1)

ax1.set_ylim(0.2, -0.2)

# Plot PACF 
plot_pacf(data_train['cnt'],lags=24, zero=False, ax=ax2)

ax2.set_ylim(0.2, -0.2)
plt.show()


result = seasonal_decompose(df_time['cnt'], model='additive', period=24*7)

# Maken van de plot
plt.figure(figsize=(15, 5))
result.seasonal.plot()
plt.title('Seizoens Component')
plt.xlim('2012-04-01', '2012-05-01')
plt.xlabel('Datum')
plt.ylabel('Seizoens-waarde')
plt.show()


def data_voorbereiding(base_df, dt_kol):
    """
    Een functie om een csv in te lezen en verschillende
    toepassingen uit te voeren op aanwezige timeseries
    data.

    Parameters:
    ----------
    base_df : pd.DataFrame
        De naam van het dataframe waar de extra
        timeseries kolommen worden gemaakt.
    
    dt_kol : str
        De naam van de kolom met de Timeseries data.

    Returns:
    ----------
    df : pd.DataFrame
        Een pandas.DataFrame waarin verschillende extra
        kolommen zijn toegevoegd waarin timeseries is
        toegepast.
    """
    
    # Omzetten van date_hour naar datetime
    df = base_df.copy()
    df[f'{dt_kol}'] = pd.to_datetime(df[f'{dt_kol}'])

    # Toegevoegd op basis van eigen ingeving
    df['Jaar'] = df[f'{dt_kol}'].dt.year.astype(int)
    df['Maand'] = df[f'{dt_kol}'].dt.month.astype(int)
    df['Uur'] = df[f'{dt_kol}'].dt.hour.astype(int)

    # Toegevoegd na aanleiding FFT
    df['is_weekend'] = (df[f'{dt_kol}'].dt.dayofweek >= 5).astype(int)

    # Aanmaken fourier voor W en D
    fourier_w = CalendarFourier(freq="W", order=1)
    fourier_d = CalendarFourier(freq="D", order=1)

    # Index datetime maken
    base_df[f'{dt_kol}'] = pd.to_datetime(base_df[f'{dt_kol}'])
    base_df = base_df.set_index(f'{dt_kol}')

    # Uitvoeren van DeterministicProcess
    dp = DeterministicProcess(
        index=base_df.index,
        constant=False,
        order=1,
        seasonal=False,
        additional_terms=[fourier_w, fourier_d],
        drop=True,
    )

    # Aanmaken van fourier dataframe
    fourier_feat = dp.in_sample()

    # Mergen van dataframes voor alle data
    df_fou = pd.merge(
        df, fourier_feat[['sin(1,freq=D)', 'cos(1,freq=D)',
                          'sin(1,freq=W-SUN)', 'cos(1,freq=W-SUN)']],
        left_on='date_hour', right_index=True)

    return df_fou


# Toepassen data_voorbereiding
df = data_voorbereiding(data_train, 'date_hour')
df.head()


# Opnieuw aanmaken Grafieken class
graf = Grafieken(df, 'cnt')

# Toepassen timeplot method
graf.timeplot('Jaar', 'weathersit')


# Toepassen timeplot method
graf.timeplot('Uur', 'is_weekend')


df_test = data_voorbereiding(data_test, 'date_hour')
df_test = df_test.set_index('date_hour').asfreq('H')
df_test


def submissie(model, suffix):
    """
    Deze functie maakt de voorspelling op de test dataset
    en vormt deze om tot een csv bestand om in te kunnen
    leveren op Kaggle.

    Parameters:
    ----------
    model : Machine Learning model
        De naam die is gegeven aan het ML-model dat
        wordt gebruikt om te voorspellen

    suffix : str
        De laatste 'tag' voor de naam van het csv
        bestand, zodat deze makkelijk te identificeren
        is na de submission.

    Returns:
    ----------
    None
        In plaats van een return maakt het een bestand
        aan in de map Kaggle Submissions. De namen zien
        er als volgt uit:

    """
    # Fitten, voorspellen en veranderen naar integer
    y_pred = model.fit(X, y).predict(df_test)
    y_pred = y_pred.astype(int)

    # Aanmaken df met alleen p_id en Outcome
    test_predictions_df = pd.DataFrame(
        {'date_hour': data_test['date_hour'],
        'cnt': y_pred})
    
    # Aanmaken van tijd
    tijd = datetime.datetime.now().strftime("%m%d%H%M%S")

    # Aanmaken csv bestand met timestamp
    test_predictions_df.to_csv(
        f'Kaggle Submissions/vs_{suffix}_{tijd}.csv',
        index=False)
    
    # print voor conformatie
    print(f'vs_{suffix}_{tijd}.csv has been saved!')

def model_score(model):
    """
    Een functie die de scores voor een model
    berekend en deze toont.

    Parameters:
    ----------
    model : elk ML model
    De naam van het gemaakte ML model

    Returns:
    ----------
    f1 : int
    De f1 score van het gegeven model
    """

    # Fitten en voorspellen met model
    y_pred = model.fit(X_train, y_train).predict(X_test)
    rmse = np.sqrt(mean_squared_error(y_test, y_pred))

    # Tonen van de resultaten
    return print(f'Root Mean Squared Error: {rmse}')

def grid_score(estimator, param_grid):
    """
    Een functie die de scores voor een grid search
    berekend en deze toont. Hierbij worden ook de
    optimale parameters getoont.

    Parameters:
    ----------
    estimator : ML model
    Het model waarop de grid search wordt toegepast

    param_grid : dict
    Een dictionary waarbij de keys parameter namen
    zijn van het model en de values verschillende
    waarden die de parameter aan kan nemen.

    Returns:
    ----------
    None
    Returned twee prints met de optimale parameters
    en beste f1 score
    """

    # Uitvoeren van de grid search
    gs = GridSearchCV(estimator=estimator,
                    param_grid=param_grid,
                    cv=cv,
                    scoring='neg_root_mean_squared_error',
                    n_jobs=-1)
    
    # Fitten van de grid search
    gs.fit(X_train, y_train)

    # Tonen van de beste score en parameters
    print(f"Beste RMSE: {abs(gs.best_score_)}")
    print(f"Beste parameters:\n{gs.best_params_}")

def feature_scan(model):
    """
    Een functie die een grafie met de feature importances kan maken.
    Dit werkt in ieder geval voor de RF en XGB modellen.

    Parameters:
    ----------
    model : ML-model
    Het model waarvan je de feature importances wilt tonen.

    Returns:
    ----------
    None
    Laat een grafiek zien met daarin de invloed van alle features.
    """
    # Aanmaken en sorteren invloed kolommen
    invloed = pd.Series(model.feature_importances_, index=X.columns)
    invloed = invloed.sort_values(ascending=True)

    # Plotten grafiek met invloed
    invloed.plot(kind='barh', figsize=(10, 6))
    plt.ylabel('Features')
    plt.xlabel('Invloed')
    plt.show()


# Aanmaken SEED voor random_state
SEED = 42

# Aanmaken van X en y
df_train = df.set_index('date_hour').asfreq('H')
X = df_train.drop('cnt', axis=1)
y = df_train['cnt']

# Splitten van de data in train en test
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=456, random_state=SEED)

# Aanmaken van de folds
cv = TimeSeriesSplit(n_splits=5, test_size=456)


# Opstellen van Linear Regression model
lr = LinearRegression(n_jobs=-1)

# Bekijken van de RMSE
model_score(lr)

Root Mean Squared Error: 117.2269334691733


param_lr = {
    'fit_intercept': [True, False],
    'copy_X' : [True, False],
    'positive' : [False, True]
}

 
# Uitvoeren van Grid Search CV
grid_score(lr, param_lr)

Beste RMSE: 126.75350985874736
Beste parameters:
{'copy_X': True, 'fit_intercept': True, 'positive': False}


# Toepassen van optimale parameters
lr_tuned = LinearRegression(n_jobs=-1)

# Bekijken van de RMSE
model_score(lr_tuned)

Root Mean Squared Error: 117.2269334691733


#submissie(lr_tuned, 'LR')


# opstellen K-Nearest Neighbors model
knn = KNeighborsRegressor(n_jobs=-1)

model_score(knn)

Root Mean Squared Error: 53.38359939354527


# Bekijken aantal neighbors
train_accuracies = {}
test_accuracies = {}
neighbors = range(1, 11)
for neighbor in neighbors:
    knn = KNeighborsRegressor(n_neighbors=neighbor)
    knn.fit(X_train, y_train)
    train_accuracies[neighbor] = knn.score(X_train, y_train)
    test_accuracies[neighbor] = knn.score(X_test, y_test)
display(train_accuracies, test_accuracies)

{1: 0.9999815166962486,
 2: 0.9643702994604253,
 3: 0.95209467887222,
 4: 0.947853355542497,
 5: 0.9451850853115388,
 6: 0.9418743474241642,
 7: 0.9390014376078347,
 8: 0.9362817072430636,
 9: 0.9333881860709422,
 10: 0.9310182177427528}

{1: 0.8355177220420447,
 2: 0.8832341262210887,
 3: 0.8887120763203489,
 4: 0.8936747310379309,
 5: 0.897679396895329,
 6: 0.9028353960226173,
 7: 0.9100914424354014,
 8: 0.9107787975884452,
 9: 0.907995193454967,
 10: 0.9083116354148302}


knn = KNeighborsRegressor(n_neighbors=4)

model_score(knn)

Root Mean Squared Error: 54.4182475724793


# parameters voor de GridSearch
params_knn = {'weights':['uniform', 'distance'],
             'metric':['minkowski','euclidean','manhattan' ],
             'algorithm':['ball_tree', 'kd_tree', 'brute']}

grid_score(knn, params_knn)

Beste RMSE: 55.75803076569273
Beste parameters:
{'algorithm': 'ball_tree', 'metric': 'minkowski', 'weights': 'uniform'}


# uitkomsten GridSearch toepassen
knn_tuned = KNeighborsRegressor(algorithm='ball_tree',
                                n_jobs=-1,
                                metric='minkowski',
                                n_neighbors=4,
                                weights='uniform')

model_score(knn_tuned)

Root Mean Squared Error: 54.4182475724793


# Opstellen van Dicision Tree Regressor model
dt = DecisionTreeRegressor(random_state=42)  

# Bekijken van de RMSE
model_score(dt)

Root Mean Squared Error: 54.01557199126803


param_dtr = {
    'criterion': ['squared_error', 'friedman_mse'],
    'max_depth': [None, 10, 20, 30],
    'min_samples_split': [2, 5],
    'min_samples_leaf': [1, 2, 4],
}

# Uitvoeren van Grid Search CV
grid_score(dt, param_dtr)

Beste RMSE: 56.58902840386772
Beste parameters:
{'criterion': 'squared_error', 'max_depth': 20, 'min_samples_leaf': 4, 'min_samples_split': 2}


# Toepassen van optimale parameters
dt_tuned = DecisionTreeRegressor(criterion='squared_error',
                                 max_depth=None,
                                 min_samples_leaf=4,
                                 min_samples_split=2)

# Bekijken van de RMSE
model_score(dt_tuned)

Root Mean Squared Error: 49.20313045211243


feature_scan(dt_tuned)


#submissie(dt_tuned, 'DT')


# Opstellen van Random Forest Regressor model
rf = RandomForestRegressor(random_state=SEED,
                              n_jobs=-1)

# Bekijken van de RMSE
model_score(rf)

Root Mean Squared Error: 37.84325840402451


param_rf = {
    'n_estimators' : [50, 100, 200],
    'max_depth' : [10, 20, 30],
    'min_samples_leaf' : [1, 2, 4]
}

grid_score(rf, param_rf)

Beste RMSE: 45.310138190556685
Beste parameters:
{'max_depth': 30, 'min_samples_leaf': 1, 'n_estimators': 200}


# Toepassen van optimale parameters
rf_tuned = RandomForestRegressor(n_estimators=200,
                           max_depth=30,
                           min_samples_leaf=1,
                           random_state=SEED,
                           n_jobs=-1)

# Bekijken van de RMSE
model_score(rf_tuned)

Root Mean Squared Error: 37.78426103437784


feature_scan(rf_tuned)


#submissie(rf_tuned, 'RF_2fa')


# Opstellen van XGBoost Regressor model
xgb = XGBRegressor(n_jobs=-1,
                   random_state=SEED)

# Bekijken van de RMSE
model_score(xgb)

Root Mean Squared Error: 37.458258453583944


# Opstellen van de parameter dict
param_xgb = {
    'n_estimators': [100, 200, 300],
    'learning_rate': [0.01, 0.1, 1.0],
    'max_depth': [5, 10, 20]
}

# Uitvoeren van de gridsearch
grid_score(xgb, param_xgb)

Beste RMSE: 41.866621142492704
Beste parameters:
{'learning_rate': 0.1, 'max_depth': 5, 'n_estimators': 300}


# Toepassen van optimale parameters
xgb_tuned = XGBRegressor(n_estimators=300,
                         learning_rate=0.1,
                         max_depth=5,
                         random_state=SEED,
                         n_jobs=-1)

# Bekijken van de RMSE
model_score(xgb_tuned)

Root Mean Squared Error: 36.402051993205


feature_scan(xgb_tuned)


# submissie(xgb_tuned, 'XGB_new')


# Opzetten basis- en meta-modellen
basis_modellen = [
    ('knn', knn_tuned),
    ('rf', rf_tuned),
    ('lr', lr_tuned)
]

meta_model = xgb_tuned

# Toepassen van StackingRegressor
sr = StackingRegressor(estimators=basis_modellen,
                       final_estimator=meta_model,
                       n_jobs=-1)

# Bekijken van de RMSE
model_score(sr)

Root Mean Squared Error: 40.11728088527348


# Opstellen van parameter dict
param_sr = {
    'passthrough' : [True, False],
    'cv' : [5, 10, None]
}

# Uitvoeren van gridsearch
grid_score(sr, param_sr)

Beste RMSE: 42.15045747739301
Beste parameters:
{'cv': 5, 'passthrough': True}


# Opzetten basis- en meta-modellen
basis_modellen = [
    ('knn', knn_tuned),
    ('rf', rf_tuned),
    ('lr', lr_tuned)
]

meta_model = xgb_tuned

# Toepassen van StackingRegressor
sr_tuned = StackingRegressor(estimators=basis_modellen,
                       final_estimator=meta_model,
                       cv=5,
                       passthrough=True,
                       n_jobs=-1)

# Bekijken van de RMSE
model_score(sr_tuned)

Root Mean Squared Error: 37.7191700700007


# submissie(sr_tuned, 'SR')


# Eerste model: Lineaire Regressie
lr_hm = LinearRegression(n_jobs=-1).fit(X_train, y_train)

# Voorspellen van de waarden voor train en test
y_pred_lr = lr_hm.predict(X_train)
pred_test = lr_hm.predict(X_test)


# Berekenen van de residuals
residuals = y_train - y_pred_lr


# Tweede model: XGBoost Regressor
xgb_hm = XGBRegressor(n_estimators=300,
                         learning_rate=0.1,
                         max_depth=5,
                         random_state=SEED,
                         n_jobs=-1).fit(X_train, residuals)

# Voorspellen van de residuals
y_pred_xgb = xgb_hm.predict(X_test)


# Samenvoegen van voorspellingen
pred = np.round(pred_test + y_pred_xgb).astype(int)

# Evalueren van hybride model
mse = mean_squared_error(y_test, pred)
print(f'Root Mean Squared Error: {np.sqrt(mse)}')

Root Mean Squared Error: 36.66069528728783


feature_scan(xgb_hm)


# Voorspellen op df_test
lr_test = lr_hm.predict(df_test)
xgb_test = xgb_hm.predict(df_test)

# Samenvoegen voorspelling
y_pred = np.round(lr_test + xgb_test).astype(int)

# Aanmaken df met alleen p_id en Outcome
test_predictions_df = pd.DataFrame(
        {'date_hour': data_test['date_hour'],
        'cnt': y_pred})

# Aanmaken van tijd
# tijd = datetime.datetime.now().strftime("%m%d%H%M%S")

# Aanmaken csv bestand met timestamp
# test_predictions_df.to_csv(
#         f'Kaggle Submissions/vs_HML_try_{tijd}.csv',
#         index=False)

# print voor conformatie
# print(f'vs_HML_try_{tijd}.csv has been saved!')


#Augmented Dickey–Fuller test:
print('Results of Dickey Fuller Test:')
dftest = adfuller(df_train['cnt'], autolag='AIC')

dfoutput = pd.Series(dftest[0:4], index=['Test Statistic','p-value','#Lags Used','Number of Observations Used'])
for key,value in dftest[4].items():
    dfoutput['Critical Value (%s)'%key] = value
    
print(dfoutput)

Results of Dickey Fuller Test:
Test Statistic                -6.625165e+00
p-value                        5.901289e-09
#Lags Used                     4.300000e+01
Number of Observations Used    1.659300e+04
Critical Value (1%)           -3.430744e+00
Critical Value (5%)           -2.861714e+00
Critical Value (10%)          -2.566863e+00
dtype: float64


fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12, 8))

# Seasonal ACF plot
plot_acf(df_time['cnt'], lags=40, ax=ax1)
ax1.set_title('Seasonal Autocorrelation Function (ACF)')
ax1.set_ylim(0.2, -0.2)

# Seasonal PACF plot
plot_pacf(df_time['cnt'], lags=40, ax=ax2)
ax2.set_title('Seasonal Partial Autocorrelation Function (PACF)')
ax2.set_ylim(0.2, -0.2)

plt.show()


#Het splitsen van de data. Dit gaat anders dan bij het splitsen van niet-tijdreeksgegevens.
split_point = len(df_train) - (3 * 7 * 24)
train, test = df_train[0:split_point], df_train[split_point:]


#Het model trainen
sarimax_model = SARIMAX(train['cnt'], order=(5, 1, 6), seasonal_order=(0, 1, 1, 24))
sarimax_model = sarimax_model.fit(disp=True, maxiter=1000)

c:\Users\crazy\AppData\Local\Programs\Python\Python312\Lib\site-packages\statsmodels\base\model.py:607: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals
  warnings.warn("Maximum Likelihood optimization failed to "


sarimax_model.plot_diagnostics(figsize=(15,12))
plt.show()


forecast = sarimax_model.get_forecast(steps=len(test))
predictions = forecast.predicted_mean
confidence_intervals = forecast.conf_int()
lower_limits = confidence_intervals.iloc[:, 0]
upper_limits = confidence_intervals.iloc[:, 1]

# Plotting
plt.figure(figsize=(12, 6))
plt.plot(train.index, train['cnt'], label='Train')
plt.plot(test.index, test['cnt'], label='Test', color = 'g')
plt.plot(predictions.index, predictions, label='Predictions', color='r')

# Fill area between confidence intervals
plt.fill_between(predictions.index, lower_limits, upper_limits, color='pink', alpha=0.5)

plt.legend()
plt.title('SARIMA Predictions with 95% Confidence Intervals')
plt.ylim(-50, 900)
plt.xlim(pd.to_datetime('2012-11-07'), pd.to_datetime('2012-11-30'))
plt.show()


# Assuming you have a test dataset called 'test' with the corresponding 'cnt' values
test['predicted'] = sarimax_model.predict(start=test.index[0], end=test.index[-1], dynamic=True)

# Calculate the squared differences between predicted and actual values
squared_diff = (test['cnt'] - test['predicted'])**2

# Calculate the mean squared error
mse = np.mean(squared_diff)

# Calculate the root mean squared error
rmse = np.sqrt(mse)

print(f'Root Mean Squared Error (RMSE): {rmse}')

Root Mean Squared Error (RMSE): 128.71632798640323

C:\Users\crazy\AppData\Local\Temp\ipykernel_23884\2082032065.py:2: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  test['predicted'] = sarimax_model.predict(start=test.index[0], end=test.index[-1], dynamic=True)


# Model trainen op gehele train set
sarimax = SARIMAX(df_train['cnt'], order=(5, 1, 6), seasonal_order=(0, 1, 1, 24))
sarimax = sarimax.fit(disp=True, maxiter=1000)


# Haal de voorspellingen op
sar_pred = sarimax_model.get_forecast(steps=len(df_test))
sar_pred = sar_pred.predicted_mean

# Maak de voorspellingen tot lijst
sar_pred = sar_pred.tolist()

# Aanmaken df voor csv
test_predictions_df = pd.DataFrame(
        {'date_hour': data_test['date_hour'],
        'cnt': sar_pred})

# Aanmaken csv
test_predictions_df.to_csv(
        f'Kaggle Submissions/vs_SARIMAX.csv',
        index=False)


#kolom namen aanpasse naar datastamp (ds) en y (target value)
df1 = df.rename(columns={'date_hour': 'ds', 'cnt': 'y'})

# Opstellen van Prophet model
prophet_model = Prophet()

# model fitten
prophet_model.fit(df1)

# dataframe maken voor predictions
# periods is 24 want seasonal period is 24 uur
future = prophet_model.make_future_dataframe(periods=24)  

# voorspelling maken
forecast = prophet_model.predict(future)

# voorspelling plotten
fig = prophet_model.plot(forecast)
plt.title('Prophet Forecast')
plt.show()

21:52:06 - cmdstanpy - INFO - Chain [1] start processing
21:52:08 - cmdstanpy - INFO - Chain [1] done processing
c:\Users\crazy\AppData\Local\Programs\Python\Python312\Lib\site-packages\prophet\plot.py:72: FutureWarning: The behavior of DatetimeProperties.to_pydatetime is deprecated, in a future version this will return a Series containing python datetime objects instead of an ndarray. To retain the old behavior, call `np.array` on the result
  fcst_t = fcst['ds'].dt.to_pydatetime()
c:\Users\crazy\AppData\Local\Programs\Python\Python312\Lib\site-packages\prophet\plot.py:73: FutureWarning: The behavior of DatetimeProperties.to_pydatetime is deprecated, in a future version this will return a Series containing python datetime objects instead of an ndarray. To retain the old behavior, call `np.array` on the result
  ax.plot(m.history['ds'].dt.to_pydatetime(), m.history['y'], 'k.',


fig2 = prophet_model.plot_components(forecast)
plt.show()

c:\Users\crazy\AppData\Local\Programs\Python\Python312\Lib\site-packages\prophet\plot.py:228: FutureWarning: The behavior of DatetimeProperties.to_pydatetime is deprecated, in a future version this will return a Series containing python datetime objects instead of an ndarray. To retain the old behavior, call `np.array` on the result
  fcst_t = fcst['ds'].dt.to_pydatetime()
c:\Users\crazy\AppData\Local\Programs\Python\Python312\Lib\site-packages\prophet\plot.py:397: FutureWarning: The behavior of DatetimeProperties.to_pydatetime is deprecated, in a future version this will return a Series containing python datetime objects instead of an ndarray. To retain the old behavior, call `np.array` on the result
  artists += ax.plot(df_y['ds'].dt.to_pydatetime(), seas[name], ls='-',
c:\Users\crazy\AppData\Local\Programs\Python\Python312\Lib\site-packages\prophet\plot.py:401: FutureWarning: The behavior of DatetimeProperties.to_pydatetime is deprecated, in a future version this will return a Series containing python datetime objects instead of an ndarray. To retain the old behavior, call `np.array` on the result
  df_y['ds'].dt.to_pydatetime(), seas[name + '_lower'],
c:\Users\crazy\AppData\Local\Programs\Python\Python312\Lib\site-packages\prophet\plot.py:397: FutureWarning: The behavior of DatetimeProperties.to_pydatetime is deprecated, in a future version this will return a Series containing python datetime objects instead of an ndarray. To retain the old behavior, call `np.array` on the result
  artists += ax.plot(df_y['ds'].dt.to_pydatetime(), seas[name], ls='-',
c:\Users\crazy\AppData\Local\Programs\Python\Python312\Lib\site-packages\prophet\plot.py:401: FutureWarning: The behavior of DatetimeProperties.to_pydatetime is deprecated, in a future version this will return a Series containing python datetime objects instead of an ndarray. To retain the old behavior, call `np.array` on the result
  df_y['ds'].dt.to_pydatetime(), seas[name + '_lower'],


actual = df['cnt']

# 24 * 365 = 8760

predicted = forecast['yhat'].tail(8760).reset_index(drop=True)

# Bekijken van de RMSE
rmse = np.sqrt(mean_squared_error(actual[-8760:], predicted)) 
print(f'RMSE: {rmse}')

RMSE: 132.69287924706384


# Train dataset
df_train_reset = df_train.reset_index()
train_data = df_train_reset.rename(columns={'date_hour': 'ds', 'cnt': 'y'})

# Maken model
proph = Prophet()
proph.fit(train_data)

# Test dataset
test_data = df_test.reset_index()
prediction_data = test_data.rename(columns={'date_hour': 'ds'})

# Voorspellingen maken
forecast = proph.predict(prediction_data)

# Aanmaken df voor csv
test_predictions_df = pd.DataFrame(
        {'date_hour': data_test['date_hour'],
        'cnt': forecast['yhat']})

# Aanmaken csv
test_predictions_df.to_csv(
        f'Kaggle Submissions/vs_Prophet.csv',
        index=False)

22:13:19 - cmdstanpy - INFO - Chain [1] start processing
22:13:22 - cmdstanpy - INFO - Chain [1] done processing

Teamleden	Kaggle Username	GitHub Username
Busse Heemskerk	bussejheemskerk	BJHeemskerk
Declan van den Hoek	declanvdh	DeclanvandenHoek
Isa Dijkstra	isadijkstra	IsaD01

	date_hour	weathersit	temp	atemp	hum	cnt
0	2011-01-01 00:00:00	1	0.24	0.2879	0.81	16
1	2011-01-01 01:00:00	1	0.22	0.2727	0.80	40
2	2011-01-01 02:00:00	1	0.22	0.2727	0.80	32
3	2011-01-01 03:00:00	1	0.24	0.2879	0.75	13
4	2011-01-01 04:00:00	1	0.24	0.2879	0.75	1

	date_hour	holiday	weathersit	temp	atemp	hum	windspeed	cnt
count	16637	16637.000000	16637.000000	16637.000000	16637.000000	16637.000000	16637.000000	16637.000000
mean	2011-12-18 02:20:08.294764544	0.028671	1.415580	0.504745	0.482608	0.624756	0.190310	190.477009
min	2011-01-01 00:00:00	0.000000	1.000000	0.020000	0.000000	0.000000	0.000000	1.000000
25%	2011-06-27 05:00:00	0.000000	1.000000	0.340000	0.333300	0.470000	0.104500	41.000000
50%	2011-12-18 06:00:00	0.000000	1.000000	0.520000	0.500000	0.620000	0.194000	143.000000
75%	2012-06-09 02:00:00	0.000000	2.000000	0.660000	0.621200	0.780000	0.253700	282.000000
max	2012-11-30 23:00:00	1.000000	4.000000	1.000000	1.000000	1.000000	0.850700	977.000000
std	NaN	0.166885	0.637298	0.192369	0.171557	0.193227	0.121915	182.026755

	date_hour	holiday	weathersit	temp	atemp	hum	windspeed	cnt
count	16637	16637.000000	16637.000000	16637.000000	16637.000000	16637.000000	16637.000000	16637.000000
mean	2011-12-18 02:20:08.294764544	0.028671	1.415580	0.504745	0.482608	0.624756	0.223709	190.477009
min	2011-01-01 00:00:00	0.000000	1.000000	0.020000	0.000000	0.000000	0.000000	1.000000
25%	2011-06-27 05:00:00	0.000000	1.000000	0.340000	0.333300	0.470000	0.122840	41.000000
50%	2011-12-18 06:00:00	0.000000	1.000000	0.520000	0.500000	0.620000	0.228047	143.000000
75%	2012-06-09 02:00:00	0.000000	2.000000	0.660000	0.621200	0.780000	0.298225	282.000000
max	2012-11-30 23:00:00	1.000000	4.000000	1.000000	1.000000	1.000000	1.000000	977.000000
std	NaN	0.166885	0.637298	0.192369	0.171557	0.193227	0.143312	182.026755

	date_hour	holiday	weathersit	temp	atemp	hum	windspeed	cnt
178	2011-01-08 17:00:00	0	1	0.16	0.1212	0.37	0.649112	69
194	2011-01-09 09:00:00	0	1	0.12	0.0758	0.46	0.614083	19
196	2011-01-09 11:00:00	0	1	0.16	0.1212	0.40	0.614083	49
265	2011-01-12 12:00:00	0	1	0.20	0.1515	0.47	0.684260	55
271	2011-01-12 18:00:00	0	1	0.20	0.1515	0.47	0.614083	137
...	...	...	...	...	...	...	...	...
15876	2012-10-28 17:00:00	0	3	0.46	0.4545	0.77	0.614083	225
16207	2012-11-13 01:00:00	0	3	0.44	0.4394	0.88	0.754438	5
16468	2012-11-23 22:00:00	0	1	0.44	0.4394	0.33	0.614083	79
16472	2012-11-24 02:00:00	0	1	0.32	0.2727	0.39	0.719290	32
16482	2012-11-24 12:00:00	0	2	0.30	0.2576	0.36	0.719290	199

Voorspellen van Verhuuraantallen¶

Inhoudsopgave ¶

H1: Inladen libaries en data ¶

H2: Explorative Data Analysis¶

H2.1: Data Opschonen ¶

H2.2: Grafieken zonder Time Series ¶

H2.3: Time Series analyses ¶

H2.3.1: Seizoens Analyse (met FFT) ¶

H2.3.2: Trend Analyse ¶

H2.3.3: Autocorrelatie en Cycli ¶

H2.4: Grafieken met Time Series ¶

H3: Modelleren ¶

H3.1: Lineaire Regressie ¶

De werking van lineaire regressie¶

Loss functie¶

Regularisatie¶

H3.2: KNN Regressor¶

H3.3: Decision tree regressor ¶

H3.4: Random Forest Regressor ¶

H3.5: XGBoost Regressor ¶

H3.6: Ensemble, StackingRegressor ¶

H3.7: Hybride model ¶

H3.8: SARIMAX Model ¶

H3.9 Phrophet ¶

H4: Bevindingen ¶

H5: Conclusie ¶

Literatuurlijst ¶

	holiday	weathersit	temp	atemp	hum	windspeed	Jaar	Maand	Uur	is_weekend	sin(1,freq=D)	cos(1,freq=D)	sin(1,freq=W-SUN)	cos(1,freq=W-SUN)
date_hour
2012-12-01 00:00:00	0	1	0.26	0.3030	0.81	0.0000	2012	12	0	1	0.000000	1.000000	-0.974928	-0.222521
2012-12-01 01:00:00	0	1	0.26	0.3030	0.81	0.0000	2012	12	1	1	0.258819	0.965926	-0.982566	-0.185912
2012-12-01 02:00:00	0	2	0.26	0.3030	0.81	0.0000	2012	12	2	1	0.500000	0.866025	-0.988831	-0.149042
2012-12-01 03:00:00	0	2	0.26	0.2727	0.81	0.1343	2012	12	3	1	0.707107	0.707107	-0.993712	-0.111964
2012-12-01 04:00:00	0	1	0.26	0.2879	0.81	0.0896	2012	12	4	1	0.866025	0.500000	-0.997204	-0.074730
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
2012-12-19 19:00:00	0	1	0.38	0.3939	0.50	0.3881	2012	12	19	0	-0.965926	0.258819	0.593820	-0.804598
2012-12-19 20:00:00	0	1	0.36	0.3485	0.57	0.2239	2012	12	20	0	-0.866025	0.500000	0.563320	-0.826239
2012-12-19 21:00:00	0	1	0.34	0.3182	0.61	0.2239	2012	12	21	0	-0.707107	0.707107	0.532032	-0.846724
2012-12-19 22:00:00	0	1	0.34	0.3485	0.61	0.0896	2012	12	22	0	-0.500000	0.866025	0.500000	-0.866025
2012-12-19 23:00:00	0	1	0.32	0.3333	0.66	0.1343	2012	12	23	0	-0.258819	0.965926	0.467269	-0.884115

Modellen	RMSE Score
Linear Regression:	117.2269334691733
Linear Regression Tuned:	117.2269334691733
KNN:	53.38359939354527
KNN Tuned:	54.4182475724793
Decision Tree:	54.01557199126803
Decision Tree Tuned:	49.12543027929688
Random Forest:	37.84325840402451
Random Forest Tuned:	37.78426103437784
XGBoost:	37.458258453583944
XGBoost Tuned:	36.402051993205
Ensemble, Stacking:	40.11728088527348
Ensemble, Stacking Tuned:	37.7191700700007
Hybrid Model:	36.66069528728783
SARIMA(X):	128.71632798640323
Prophet:	132.69287924706384